
R version 4.1.2 (2021-11-01) -- "Bird Hippie"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: aarch64-apple-darwin20 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.


> library(support.CEs)
Registered S3 method overwritten by 'DoE.base':
  method           from       
  factorize.factor conf.design
> library(tidyverse)
── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.6     ✔ purrr   0.3.4
✔ tibble  3.1.7     ✔ dplyr   1.0.9
✔ tidyr   1.2.0     ✔ stringr 1.4.0
✔ readr   2.1.2     ✔ forcats 0.5.1
── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
> library(foreign)
> library(survival)
> library(modelsummary)
> library(cjoint)
Loading required package: sandwich
Loading required package: lmtest
Loading required package: zoo

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

Loading required package: survey
Loading required package: grid
Loading required package: Matrix

Attaching package: ‘Matrix’

The following objects are masked from ‘package:tidyr’:

    expand, pack, unpack


Attaching package: ‘survey’

The following object is masked from ‘package:graphics’:

    dotchart

cjoint: AMCE Estimator for Conjoint Experiments
Version: 2.1.0
Authors: Soubhik Barari, Elissa Berwick, Jens Hainmueller, Daniel Hopkins, Sean Liu, Anton Strezhnev, Teppei Yamamoto


Attaching package: ‘cjoint’

The following object is masked from ‘package:tibble’:

    view

> library(ggplot2)
> library(data.table)
data.table 1.14.2 using 1 threads (see ?getDTthreads).  Latest news: r-datatable.com
**********
This installation of data.table has not detected OpenMP support. It should still work but in single-threaded mode.
This is a Mac. Please read https://mac.r-project.org/openmp/. Please engage with Apple and ask them for support. Check r-datatable.com for updates, and our Mac instructions here: https://github.com/Rdatatable/data.table/wiki/Installation. After several years of many reports of installation problems on Mac, it's time to gingerly point out that there have been no similar problems on Windows or Linux.
**********

Attaching package: ‘data.table’

The following objects are masked from ‘package:dplyr’:

    between, first, last

The following object is masked from ‘package:purrr’:

    transpose

> library(readstata13)
> #plot theme
> plot_theme =  theme_bw()+ theme(text = element_text(size=10), 
+                                 panel.grid.major = element_blank(),
+                                 panel.grid.minor = element_blank(),
+                                 panel.grid.major.y = element_line(colour = "#CCCCCC"),
+                                 axis.line = element_line(colour = "black"),
+                                 strip.background=element_rect(colour=NA, fill="#DDDDDD"),
+                                 legend.title=element_blank(), 
+                                 legend.position="top",
+                                 axis.title=element_text(size=14,face="bold"), 
+                                 axis.text=element_text(size=12),
+                                 strip.text.x = element_text(size = 12),
+                                 legend.text=element_text(size=12), 
+                                 axis.text.x = element_text(angle = 45, hjust = 1, vjust=1)) 
> #plot features
> pd<-position_dodge(.2)
> pd_text<-position_dodge(1)
> #number of runs for bootstrap 
> runs = 2500
> 
> results_data <- read.dta13("lucid_formatted.dta")
Warning message:
In read.dta13("lucid_formatted.dta") : 
   Factor codes of type double or float detected in variables

   pid3_lean

   No labels have been assigned.
   Set option 'nonint.factors = TRUE' to assign labels anyway.

> 
> results_data<-results_data %>% 
+   select(id, profile_chosen,mark_up, country, rating, task_num, ethnocentrism) %>% 
+   rename(respondent_id = id)
> 
> #format
> results_data<- results_data %>% 
+   mutate(country = factor(country, levels=c("United States", "India", "China","Canada","Japan"))) %>%
+   mutate(mark_up_cat = case_when(
+     mark_up == "1" ~ "100", 
+     mark_up == "0.25" ~ "25", 
+     mark_up == "0.5" ~ "50",
+     mark_up == "0" ~ "0")) %>% 
+   mutate(mark_up_cat =factor(mark_up_cat, levels=c("0", "25", "50", "100"))) %>% 
+   mutate(rating_cat = factor(rating, levels =c("5 out of 5 stars", "4 out of 5 stars", "3 out of 5 stars")))
> 
> 
> with(results_data, summary(country))
United States         India         China        Canada         Japan 
         5249          5248          5294          5397          5268 
> with(results_data, summary(mark_up_cat))
   0   25   50  100 
6599 6583 6633 6641 
> 
> ##ethno_bins (new ethnocentrism measure)
> with(results_data, summary(ethnocentrism))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.4375  0.5625  0.5578  0.6875  1.0000 
> 
> results_data <- results_data %>%
+   mutate(ethno_factor = case_when(
+     ethnocentrism <= .4375 ~ "Low ethnocentrism", 
+     ethnocentrism > .4375 & ethnocentrism < .6875  ~ "Medium ethnocentrism", 
+     ethnocentrism >= .6875  ~ "High ethnocentrism", 
+   )) %>% mutate(ethno_factor = as.factor(ethno_factor))
> 
> with(results_data, table(ethno_factor))
ethno_factor
  High ethnocentrism    Low ethnocentrism Medium ethnocentrism 
                7814                 7138                11504 
> 
> #################
> # marginal means
> ################
> 
> library("cregg")

Attaching package: ‘cregg’

The following object is masked from ‘package:cjoint’:

    amce

> 
> #specify our cjoint design 
> c_design <- profile_chosen ~ mark_up_cat + rating_cat + country 
> 
> # marginal means
> mm_full <- cj(results_data, c_design, id = ~respondent_id, estimate = "mm")
Warning messages:
1: In logLik.svyglm(x) : svyglm not fitted by maximum likelihood.
2: In logLik.svyglm(x) : svyglm not fitted by maximum likelihood.
3: In logLik.svyglm(x) : svyglm not fitted by maximum likelihood.
> mm_full$BY = "Full sample"
> 
> ### marginal means  -- by ethno
> mm_ethno_bins <- cj(results_data, c_design, id = ~respondent_id, by = ~ ethno_factor, estimate = "mm")
Warning messages:
1: In logLik.svyglm(x) : svyglm not fitted by maximum likelihood.
2: In logLik.svyglm(x) : svyglm not fitted by maximum likelihood.
3: In logLik.svyglm(x) : svyglm not fitted by maximum likelihood.
4: In logLik.svyglm(x) : svyglm not fitted by maximum likelihood.
5: In logLik.svyglm(x) : svyglm not fitted by maximum likelihood.
6: In logLik.svyglm(x) : svyglm not fitted by maximum likelihood.
7: In logLik.svyglm(x) : svyglm not fitted by maximum likelihood.
8: In logLik.svyglm(x) : svyglm not fitted by maximum likelihood.
9: In logLik.svyglm(x) : svyglm not fitted by maximum likelihood.
> 
> mm_full_plot  <- mm_full %>% filter(feature == "country") %>% 
+   select(BY, estimate, lower, upper, statistic, level ) 
> mm_ethno_bins_plot  <- mm_ethno_bins %>% filter(feature == "country" & BY != "Medium ethnocentrism") %>% 
+   select(BY, estimate, lower, upper, statistic, level ) 
> 
> 
> data_for_cjoint <- results_data %>%  
+   drop_na() %>%
+   mutate(country = relevel(country, ref="United States")) %>%
+   mutate(rating_cat = relevel(rating_cat, ref="5 out of 5 stars"))
> summary(data_for_cjoint)
 respondent_id    profile_chosen    mark_up                country        rating             task_num      ethnocentrism    mark_up_cat
 Min.   :   1.0   Min.   :0.0    Min.   :0.0000   United States:5249   Length:26456       Min.   : 1.000   Min.   :0.0000   0  :6599   
 1st Qu.: 281.0   1st Qu.:0.0    1st Qu.:0.2500   India        :5248   Class :character   1st Qu.: 3.000   1st Qu.:0.4375   25 :6583   
 Median : 559.0   Median :0.5    Median :0.5000   China        :5294   Mode  :character   Median : 6.000   Median :0.5625   50 :6633   
 Mean   : 559.7   Mean   :0.5    Mean   :0.4386   Canada       :5397                      Mean   : 6.478   Mean   :0.5578   100:6641   
 3rd Qu.: 840.0   3rd Qu.:1.0    3rd Qu.:1.0000   Japan        :5268                      3rd Qu.: 9.000   3rd Qu.:0.6875              
 Max.   :1117.0   Max.   :1.0    Max.   :1.0000                                           Max.   :12.000   Max.   :1.0000              
            rating_cat                 ethno_factor  
 5 out of 5 stars:8808   High ethnocentrism  : 7814  
 4 out of 5 stars:8895   Low ethnocentrism   : 7138  
 3 out of 5 stars:8753   Medium ethnocentrism:11504  
                                                     
                                                     
                                                     
> 
> with(results_data, table(ethno_factor))
ethno_factor
  High ethnocentrism    Low ethnocentrism Medium ethnocentrism 
                7814                 7138                11504 
> 
> set.seed(1234)
> results_vector = list()
> estimates_results = list()
> resp_ids <- unique(data_for_cjoint$respondent_id)
> 
> for (run in 1:runs) {
+ 
+   samp <- sample(x = resp_ids, size = length(resp_ids), replace = TRUE)
+   
+   bootSample <- as.data.table(data_for_cjoint)
+   setkey(bootSample, "respondent_id")
+   bootSample <- bootSample[J(samp), allow.cartesian = TRUE]
+   
+   #store complete results
+   results_vector[[run]] <- cjoint::amce(profile_chosen ~ mark_up_cat + rating_cat + country, data=bootSample, respondent.id = "respondent_id")
+   
+   #store estimates 
+   country <-  as.tibble(t(results_vector[[run]]$estimates$country), rownames="est")
+   markup <- as.tibble(t(results_vector[[run]]$estimates$markupcat), rownames="est")
+   #store the estimates
+   estimates_results[[run]] <- bind_rows(country, markup) %>% mutate(run = run)
+   
+ }
Warning message:
`as.tibble()` was deprecated in tibble 2.0.0.
Please use `as_tibble()` instead.
The signature and semantics have changed, see `?as_tibble`.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated. 
> 
> 
> full_sample_results <- bind_rows(estimates_results) 
> 
> full_sample_markups <- full_sample_results %>% filter(grepl('markup', est)) %>% 
+   select(est, AMCE, run) %>%
+   rename(markup_est = AMCE, markup = est)
> 
> 
> full_sample_results_summary_disag <-  full_sample_results %>% 
+   filter(grepl('country', est)) %>%  
+   left_join(full_sample_markups, by=c("run")) %>%
+   mutate(markup_num = str_replace(markup, "markupcat", "")) %>%
+   mutate(mwtp = AMCE/(markup_est/as.numeric(markup_num))) %>%
+   group_by(est, markup_num) %>%
+   summarise(mean_amce = mean(AMCE), 
+             lower_amce =quantile(AMCE, probs=c(.025)), 
+             upper_amce = quantile(AMCE, probs=c(.975)), 
+             mean_mwtp = mean(mwtp), 
+             lower_mwtp =quantile(mwtp, probs=c(.025)), 
+             upper_mwtp = quantile(mwtp, probs=c(.975))) %>%
+   mutate(ethno_factor = "Full sample") %>%
+   rename(main = est)
`summarise()` has grouped output by 'est'. You can override using the `.groups` argument.
> 
> 
> amce_full_sample_results_disag <- full_sample_results_summary_disag %>% select(main, mean_amce, lower_amce, upper_amce, ethno_factor, markup_num) %>%
+   rename(mean = mean_amce, upper = upper_amce, lower=lower_amce) %>% mutate(analysis = "AMCE")
> 
> mwtp_full_sample_results_disag <- full_sample_results_summary_disag %>% select(main, mean_mwtp, lower_mwtp, upper_mwtp, ethno_factor, markup_num ) %>%
+   rename(mean = mean_mwtp, upper = upper_mwtp, lower=lower_mwtp) %>% mutate(analysis = "MWTP")
> 
> 
> 
> 
> full_sample_results_summary <-  full_sample_results %>% filter(grepl('country', est)) %>%  left_join(full_sample_markups, by=c("run")) %>%
+   mutate(markup_num = str_replace(markup, "markupcat", "")) %>%
+   mutate(mwtp = AMCE/(markup_est/as.numeric(markup_num))) %>%
+   group_by(est) %>%
+   summarise(mean_amce = mean(AMCE), 
+             lower_amce =quantile(AMCE, probs=c(.025)), 
+             upper_amce = quantile(AMCE, probs=c(.975)), 
+             mean_mwtp = mean(mwtp), 
+             lower_mwtp =quantile(mwtp, probs=c(.025)), 
+             upper_mwtp = quantile(mwtp, probs=c(.975))) %>%
+   mutate(ethno_factor = "Full sample") %>%
+   rename(main = est)
> 
> amce_full_sample_results <- full_sample_results_summary %>% select(main, mean_amce, lower_amce, upper_amce, ethno_factor) %>%
+   rename(mean = mean_amce, upper = upper_amce, lower=lower_amce) %>% mutate(analysis = "AMCE")
> 
> mwtp_full_sample_results <- full_sample_results_summary %>% select(main, mean_mwtp, lower_mwtp, upper_mwtp, ethno_factor ) %>%
+   rename(mean = mean_mwtp, upper = upper_mwtp, lower=lower_mwtp) %>% mutate(analysis = "MWTP")
> 
> 
> 
> 
> ###non-parametric bootstrap 
> 
> results_vector = list()
> estimates_results = list()
> resp_ids <- unique(data_for_cjoint$respondent_id)
> set.seed(1234)
> 
> for (run in 1:runs) {
+ 
+   samp <- sample(x = resp_ids, size = length(resp_ids), replace = TRUE)
+   bootSample <- as.data.table(data_for_cjoint)
+   setkey(bootSample, "respondent_id")
+   bootSample <- bootSample[J(samp), allow.cartesian = TRUE]
+   
+   #store complete results
+   results_vector[[run]] <- cjoint::amce(profile_chosen ~ mark_up_cat + rating_cat + country*ethno_factor, data=bootSample,  respondent.varying = "ethno_factor", respondent.id = "respondent_id")
+   
+   #store estimates 
+   country_ethno <-  as.tibble(t(results_vector[[run]]$cond.estimates$`country:ethnofactor`), rownames="est")
+   country <-  as.tibble(t(results_vector[[run]]$cond.estimates$country), rownames="est")
+   ethno <- as.tibble(t(results_vector[[run]]$cond.estimates$ethnofactor), rownames="est")
+   markup <- as.tibble(t(results_vector[[run]]$estimates$markupcat), rownames="est")
+   estimates_results[[run]] <- bind_rows(country_ethno, country, ethno, markup) %>% mutate(run = run)
+   
+ }
> 
> 
> results <- bind_rows(estimates_results) %>% separate(est, sep=":", c("main", "conditional") )
Warning message:
Expected 2 pieces. Missing pieces filled with `NA` in 22500 rows [9, 10, 11, 12, 13, 14, 15, 16, 17, 26, 27, 28, 29, 30, 31, 32, 33, 34, 43, 44, ...]. 
> 
> markups <- results %>% filter(grepl('markup', main)) %>% 
+   select(main, AMCE, run) %>%
+   rename(markup_est = AMCE, markup = main)
> 
> high_ethno <- results %>% filter(grepl('country', main) & is.na(conditional)) %>% 
+   mutate(ethno_factor = "High ethnocentrism") %>%
+   select(`Conditional Estimate`, ethno_factor, run, main) %>%
+   rename(estimate = `Conditional Estimate` )
> 
> low_ethno <- results %>% 
+   filter(grepl('country', main) & (is.na(conditional) | conditional == "ethnofactorLowethnocentrism") ) %>%
+   mutate(ethno_factor = "Low ethnocentrism") %>%
+   group_by(run, main, ethno_factor) %>%
+   summarize( estimate =  sum (`Conditional Estimate`))  %>%
+   ungroup()
`summarise()` has grouped output by 'run', 'main'. You can override using the `.groups` argument.
>   
> 
> 
> full_results_disag <- bind_rows(high_ethno, low_ethno) %>% 
+   left_join(markups, by=c("run")) %>%
+   mutate(markup_num = str_replace(markup, "markupcat", "")) %>%
+   mutate(mwtp = estimate/(markup_est/as.numeric(markup_num))) %>%
+   group_by(main, ethno_factor, markup_num) %>%
+   summarise(mean_amce = mean(estimate), 
+             lower_amce =quantile(estimate, probs=c(.025)), 
+             upper_amce = quantile(estimate, probs=c(.975)), 
+             mean_mwtp = mean(mwtp), 
+             lower_mwtp =quantile(mwtp, probs=c(.025)), 
+             upper_mwtp = quantile(mwtp, probs=c(.975)))
`summarise()` has grouped output by 'main', 'ethno_factor'. You can override using the `.groups` argument.
> 
> 
> 
> amce_results_disag <- full_results_disag %>% select(main, mean_amce, lower_amce, upper_amce, ethno_factor, markup_num) %>%
+   rename(mean = mean_amce, upper = upper_amce, lower=lower_amce) %>% mutate(analysis = "AMCE")
> 
> mwtp_results_disag <- full_results_disag %>% select(main, mean_mwtp, lower_mwtp, upper_mwtp, ethno_factor, markup_num) %>%
+   rename(mean = mean_mwtp, upper = upper_mwtp, lower=lower_mwtp) %>% mutate(analysis = "MWTP")
> 
> 
> full_results <- bind_rows(high_ethno, low_ethno) %>% 
+   left_join(markups, by=c("run")) %>%
+   mutate(markup_num = str_replace(markup, "markupcat", "")) %>%
+   mutate(mwtp = estimate/(markup_est/as.numeric(markup_num))) %>%
+   group_by(main, ethno_factor) %>%
+   summarise(mean_amce = mean(estimate), 
+             lower_amce =quantile(estimate, probs=c(.025)), 
+             upper_amce = quantile(estimate, probs=c(.975)), 
+             mean_mwtp = mean(mwtp), 
+             lower_mwtp =quantile(mwtp, probs=c(.025)), 
+             upper_mwtp = quantile(mwtp, probs=c(.975)))
`summarise()` has grouped output by 'main'. You can override using the `.groups` argument.
> 
> 
> amce_results <- full_results %>% select(main, mean_amce, lower_amce, upper_amce, ethno_factor) %>%
+   rename(mean = mean_amce, upper = upper_amce, lower=lower_amce) %>% mutate(analysis = "AMCE")
> 
> mwtp_results <- full_results %>% select(main, mean_mwtp, lower_mwtp, upper_mwtp, ethno_factor ) %>%
+   rename(mean = mean_mwtp, upper = upper_mwtp, lower=lower_mwtp) %>% mutate(analysis = "MWTP")
> 
> mwtp_results
# A tibble: 8 × 6
# Groups:   main [4]
  main           mean lower upper ethno_factor       analysis
  <chr>         <dbl> <dbl> <dbl> <chr>              <chr>   
1 countryCanada  48.4  34.9  64.1 High ethnocentrism MWTP    
2 countryCanada  32.6  21.4  45.4 Low ethnocentrism  MWTP    
3 countryChina  104.   81.9 130.  High ethnocentrism MWTP    
4 countryChina   82.8  63.0 105.  Low ethnocentrism  MWTP    
5 countryIndia   83.4  64.0 106.  High ethnocentrism MWTP    
6 countryIndia   63.5  47.5  82.0 Low ethnocentrism  MWTP    
7 countryJapan   58.8  43.1  76.3 High ethnocentrism MWTP    
8 countryJapan   55.8  41.2  72.6 Low ethnocentrism  MWTP    
> 
> #combine the cjoint results with the cregg results
> plot_results <- bind_rows(amce_results, mwtp_results, amce_full_sample_results, mwtp_full_sample_results) %>%
+   rename(estimate = mean, BY = ethno_factor, statistic = analysis, level = main)
> 
> 
> 
> plot_results_disag <- bind_rows(amce_full_sample_results_disag, mwtp_full_sample_results_disag, amce_results_disag, mwtp_results_disag) %>%
+   rename(estimate = mean, BY = ethno_factor, statistic = analysis, level = main)
> 
> 
> 
> #analysis labels
> 
> analysis_labels <- c("Marginal mean\n(relative to all other countries)", "AMCE\n(relative to United States)", "Implied MWTP\n(relative to United States)")
> 
> for_plot <- bind_rows(plot_results, mm_full_plot, mm_ethno_bins_plot ) %>%
+   mutate(level = str_remove(level, "^country")) %>%
+   mutate(analysis_label = case_when(
+            statistic == "AMCE" ~ analysis_labels[2], 
+            statistic == "MWTP" ~analysis_labels[3],
+            statistic == "mm" ~ analysis_labels[1], 
+            TRUE ~ as.character(statistic))) %>% 
+   mutate(analysis_label = factor(analysis_label, levels = analysis_labels )) %>%
+   mutate(BY = case_when(
+     BY == "Full sample" ~ "Sample-wide average", 
+     TRUE ~ BY )) %>%
+   mutate(BY = factor(BY, levels=c("Sample-wide average", "Low ethnocentrism", "High ethnocentrism"))) %>%
+   mutate(level = factor(level, levels = c("Canada","Japan","China","India","United States"))) 
>   
>   
> 
> pd_combined<-position_dodge(.75)
> combined_plot_short_article <- ggplot(data=for_plot, aes(y=estimate, x=level, group=BY)) +
+   geom_errorbar(aes(ymin=lower, ymax=upper), width=0, position=pd_combined)+
+   geom_point(aes(fill=BY, shape=BY),color="black",size=4,position=pd_combined) + 
+   xlab("Country") +
+   ylab("Estimate") +
+   scale_shape_manual(values=c(22,21, 23, 24))+
+   facet_wrap(analysis_label~., scales="free") + plot_theme
> combined_plot_short_article
> 
> 
> ggsave("FigureA6.png", combined_plot_short_article, dpi=600, width=14, height=7)
> 